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Abstract 

We introduce a mesoscale method for simulating hydrodynamic transport and 
self assembly of inhomogeneous polymer melts in pressure driven and drag induced 
flows. This method extends dynamic self consistent field theory (DSCFT) into the 
hydrodynamic regime where bulk material transport and viscoelastic effects play 
a significant role. The method combines four distinct components as a single cou- 
pled system, including (1) non-equilibrium self consistent field theory describing 
block copolymer self-assembly, (2) multi-fluid Navier-Stokes type hydrodynamics for 
tracking material transport, (3) constitutive equations modeling viscoelastic phase 
separation, and (4) rigid wall fields which represent moving channel boundaries, ma- 
chine components, and nano-particulate fillers. We also present an efficient, pseudo- 
spectral implementation for this set of coupled equations which enables practical 
application of the model in periodic domains. We validate the model by repro- 
ducing well known phenomena including equilibrium diblock meso-phases, analytic 
Stokes flows, and viscoelastic phase separation of glassy /elastic polymer melts. We 
also demonstrate the stability and accuracy of the numerical implementation by 
examining its convergence under grid-size refinement. 



1 Introduction 

Many important products are composed of inhomogeneous polymeric fluids which 
are processed in molten form and subjected to industrial techniques such as melt 
injection, blow molding, spin casting, and fiber drawing. The physical properties of 
these products are often dictated by the detailed distribution of their constituent 
components, which in turn is determined by the manner in which they were pro- 
cessed. Improving the final product may be achieved by improving the process, and 
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computational tools can provide a convenient and cost effective alternative to trial- 
and-error refinement. At the minimum, such a numerical tool must be capable of 
simulating phase separation, polymer viscoelasticity, boundary surface wetting, the 
effects of contaminants, and the manner in which hydrodynamic transport influences 
these phenomena. 

In this paper, we propose a method for simulating the transport of inhomoge- 
neous polymeric fluids, as well as an eflicient pseudo-spectral implementation. The 
method combines a non-equilibrium version of Self- Consistent Field Theory (SCFT), 
a Navier-Stokes type hydrodynamic model, and a set viscoelastic constitutive equa- 
tions, into a single, coupled set of nonlinear differential equations. It also employs 
a set of continuous, rigid wall fields, which may be moving, to simulate pressure 
driven flows, drag induced flows, and boundary-wetting conditions. Subsequently, 
we refer to the method as Hydrodynamic Self Consistent Field Theory (HSCFT) in 
order to underscore its emphasis on the hydrodynamic transport of polymeric fluids. 
An schematic diagram illustrating the interaction of each component is presented 
in flg.l. 

In contrast to "phase field" techniques [1,2,3] which employ a Ginzburg-Landau 
free energy, this method is capable of simulating the assembly of multiblock copoly- 
mer meso-phases where the polymeric nature of the chains is explicitly taken into 
account. It also differs from traditional SCFT and DSCFT methods which are 
incapable of simulating the effects of hydrodynamic transport in pressure driven 
and drag induced flows. Generally speaking, SCFT [4,5,6] describes equilibrium 
morphologies and meso-phases boundaries, DSCFT [7,8,9,10,11,12] models non- 
equilibrium systems in which hydrodynamic transport may be neglected including 
phase separating melts and systems subjected to simple shear fields, and HSCFT is 
appropriate for the description of non-equilibrium systems in which hydrodynamic 
effects play an important role. 

In section 2 we present the governing equations for each component of the model. 
In section 3, we convert these equations to a dimcnsionless form, and in section 4 we 
present a pseudo-spectral numerical method suitable for practical implementation 
of the method on a computer. We validate the method in section 5 by reproducing 
the results of established numerical techniques and demonstating qualitative agree- 
ment with experiment and verify the stability and convergence of the numerical 
implementation under grid-size refinement. We summarize our results in section 6. 
For completeness, we also present derivations of the SCFT thermodynamic model 
in appendix A and the multi-fluid hydrodynamic model in appendix B. 

2 Governing Equations 

2.1 Polymer Thermodynamics 

The micro-physical, material specific properties of the system are modeled by 
a non-equilibrium version of SCFT. Self consistent field theory is a mean field, 
mesoscale technique ideal for simulating the thermodynamic properties of entangled 
polymers, which is capable of simulating phase separation of incompatible species, 
polymer-wall interactions, and the effects of polymer polydispersity. Following the 
formal procedure discussed in appendix A, we obtain the free energy of a blend of 
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Fig. 1. Overview of the HSCFT scheme. The thermodynamic model tracks the 
micro-physical behavior of the polymers, producing imbalanced chemical potential 
fields which, together with moving walls and external body forces, induce bulk 
hydrodynamic motion. The multi-fluid model iterates to satisfy the no flow, no 
slip, constant density, and force balance conditions, and the resultant pressure and 
velocity fields are used to transport the volume fractions. 



C distinct copolymer chains in the canonical ensemble, which takes the form 

F = U^^ + U^^ - kTS 
where the enthalpic components of the free energy are 
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The U(j)(j) term represents the energy associated with all two body monomer-monomer 
interactions, where 0i(r) is the volume fraction of monomers of type i at point r. 
The sum runs over all M distinct monomer species. Note that for our purposes we 
consider two species of the same monomer type to be distinct if they belong to differ- 
ent copolymers, as the two types will exhibit different dynamic behavior. The Flory 
incompatibility parameter, Xiji sets the strength of the effective repulsion between 
them, which is positive for immiscible polymers. Similarly, the U^^ term represents 
the net energy of all monomer- wall interactions, where i/jj{r) represents the local 
volume fraction of solid material of type j out of W possible "wall" materials. 
The free energy also contains the term 

c 1 r ^ 

S/k = - Vn«lng„[{u;}] - - / dr' Vcc;i(r')0i(r') (4) 



which describes the net entropy arising from all copolymers in the system. Each 
term Qa represents a partition function over all possible configurations of a given 
copolymer a when subjected to the externally applied fields ijJi{r) where i runs over 
all distinct monomers species and the second term "subtracts off' the effects of the 
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external field. Note that the copolymer index a is a function of monomer index 
i. For example, in a blend containing two triblock copolymer species, a = 1 for 
i G {1, 2, 3} and a = 2 for i € {4, 5, 6}. 

The manner in which the copolymers chains are stretched and arranged in a 
mean-field sense may be determined by solving two Feynman-Kac style diffusion 
equations over the polymerization index s for each copolymer a, 

dsQair, s) = R^g^V'^qair, s) - Na^a{r, s)qair, s) (5) 

dsqiir, s) = Rl^V\i{y, s) + N^^^{v, s)qi{v, s) (6) 

The function g(r, s) is called the propagator and (?^(r, s) is the co-propagator, where 
q{r , s)q\r , s) / Q represents the probability of finding segment s at the point r 
given the initial conditions (?(r, 0) = 1 and g^(r, 1) = 1. The function ria(r, s) = 
'Y^^\Oi>i{v)^i{s) is the external field acting on the polymer at s, 7i(s) is an occu- 
pation function, indicating the linear density of monomer i at index s, and Rg^ 
is the radius of gyration of polymer a in the absence of external fields. Once the 
propagator equations have been solved, the partition functions may be computed 
by averaging over the system volume V for arbitrary s' . 

Q„ = l fdr'q^{r',s')q\v',s') (7) 



Functional differentiation of the free energy with respect to each of the fields 
4>i{r), V'i(r), and LOi{r) gives the thermodynamic potentials which drive phase sep- 
aration in non-equilibrium states. 

kT I ^ W \ 
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Mr(r) = T^ = — I'd^ qair,s)qUr,sH{s) - <^,(r) ) (10) 

Vo \Qa Jo 



The chemical potentials iJ^fir) and /if (r) correspond to conserved number density 
parameters and are typically slow to relax, as transport over large distances may be 
required to achieve equilibrium. The chemical potential fields /i^, on the other hand 
correspond to the non-conserved, local rearrangement of chains, and as such may 
be expected to remain close to their equilibriTim value, /^'^(r) = 0. This condition 
implicitly defines a set of constraints on iOi{r), which may be solved numerically, in 
an iterative manner. 



2.2 Hydrodynamic Transport 

Simulating the flow of inhomogeneous polymeric fluids in industrial processing 
flows requires a hydrodynamic model that tracks multiple viscoelastic fluids and 
irregular, possibly moving channel walls. Both requirements may be satisfied simul- 
taneously by combining a generalization of the "two-fiuid model" for polymer blends 
[15] with flow penalization techniques for irregular boundary surfaces described in 
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[16]. The interested reader will find the details of this derivation in appendix B. The 
resulting multiple-fluid model is composed of a set of modified Navier-Stokes trans- 
port equations which relate the forces in the system to the transport of conserved 
quantities. 

In the absence of chemical reactions, the volume fraction fields are conserved as 
expressed by the continuity equations 

dtct^i + V ■(l>ivf = (11) 
StVi + V-Vivf = (12) 

where (r) is the velocity of fluid (?!>i(r) and v^(r) is the velocity of the solid 

material ^i(r). 

It is convenient to divide the velocity of each fluid into collective and relative 
components, vf = vt + w. The tube velocity = Hf^afvf + Sj' aj vj repre- 
sents the collective motion of the network of topological constrains in an entangled 
polymer melt, as described in Brochard's theory of mutual diffusion [28], and the 
stress division parameters a'f' and a"^ are obtained by balancing the frictional forces 
acting on the network as described in [15]. The quantity Wj = vf — vr represents 
the velocity of fluid i relative to the network. 

The relative velocity of each component may be obtained from the momentum 
transport equations in the low Reynold's number limit, 

Cf (vf - vt) = af V ■ a - ^ + Vp) + if (13) 
Cf (vf ~ vt) = atV-a- + Vp) + ff (14) 

where Cf is the friction coefficient associated with fluid i and is the friction 
coefficient associated with solid component j. Inertial effects are negligible due to the 
small size and high viscosity of the system and the fluid velocities v*^ are maintained 
in a psuedo-steady state such that frictional forces balance the osmotic forces (piVfif, 
pressure gradients Vp, viscoelastic stresses V • <t and external body forces if. If the 
velocity field of the solid objects are specified, then v^(r) is known, and Eq. (14) 
may be solved for f/'(r). Integrating this quantity over each rigidly connected region 
allows us to measure the net force and torque acting on that object. 

Summing over the relative transport equations gives the modified Navier-Stokes 
momentum transport equation in the low Reynolds number limit, 

= Vtt + Vp - V • o- - - (15) 

where we have defined the total osmotic pressure gradient Vvr = T.f(piVfif + 
T^Y ipj^ fJ'f , the net wall force f'^(r) = Sj^]^fj^(r) and the net body force f'^(r) = 
Tifl-^^if (r) . This equation may be solved simultaneously with the incompressible con- 
tinuity equation V-v = to obtain the mean velocity field v(r) = Sj^jV^' + SjV'jVj', 
and pressure field. Once v and are known, the tube velocity may obtained from 
the relationship. 

1 - 2^fe=i(% - w 

In the case where the fluids axe d3mamically matched, the stress division parameters 
reduce to the volume fractions, and the tube velocity reduces to the mean velocity. 
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2.3 Viscoelastic Constitutive Equations 

Polymers are viscoelastic materials which flow like a liquid over long times and 
behave like elastic solids when subjected to rapidly changing stresses and strains. 
Unlike simple Newtonian fluids, the clastic nature of polymers allows them to store 
energy, resulting in a time dependent stress-strain relationship which is described by 
the material's constitutive equation. While many constitutive equations are avail- 
able which approximate these viscoelastic propeties, we have chosen to employ the 
Oldyrod-B model which allows us to study a purely viscous fluid (Gj = = 0), a 
purely viscoelastic melt (rjs = 0), or anything in-between. Following Tanaka's exam- 
ple [2], we employ shear moduli of the form Gi{(p) = Goi^f and bulk moduli of the 
form Ki{(l)) = Koi6 {(pi — /,) where fi is the average volume fraction of component 
i in the system and 6 is the step function. These choices enable us to simulate sys- 
tems composed of materials with a large dynamic contrast such as a glassy elastic 
polymer blends and facilitates comparison with well established Ginzbug-Landau 
methods [2]. 

The total viscoelastic force at point r is found by summing the elastic forces con- 
tributed by each polymeric component and a dissipative viscous force as described 

by 

V-o-(r) = ^V-o-i(r) + %v2vT(r) (17) 

i 

where the shear and bulk stresses stored in component i evolve over time according 
to the constitutive equation 

= Gi{<t>)KT + Ki{(t>){V ■ vt)5 - o-^/ra (18) 

Each of the derivatives above is an upper convected derivative, a = dta + vt • Ver — 
<T ■ Vvt — (Vvt)^ • <T, which tracks the transport of elastic stresses due to fluid 
convection while eliminating spurious stresses induced by purely rotational motion, 
and the quantity kt = Vvt + (Vvt)"'" ~ | " vr) ^ is the shear strain-rate tensor. 

This constitutive equation is an appropriate for polymeric fluids subjected to 
moderate strain-rate shearing flows. For simulations with large elastic strains or 
highly extensional flows, the constitutive equation may be replaced with a more 
sophisticated phenomenological model or one based upon the SOFT microphysics 
[17]. 

3 Nondimensionalization 

We convert the governing equations to a dimensionless form by extracting a 
characteristic dimensional scale from each quantity. 

r = r/Le t = t/tc v = -v/Vc F = F/E^ a = <7/ac (19) 

The characteristic length-scale Lc is set by the radius of gyration of the smallest 

unperturbed copolymer chain, Lc = Rg = by^N/2d, where b is the statistical seg- 
ment length, d is the dimension of the system, and N is its polymerization index . 
The characteristic time-scale is dictated by the smallest reptative disentanglement 
time tc = min(T„), and the characteristic energy is chosen to be the thermal energy, 
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Ec = kT. The chartacteristic viscoelastic stress = min(Gj) is set by the smahest 
polymer shear modulus, and the characteristic convectivc velocity is the defined 
by the ratio Vc = Lc/tc- With these definitions, the propagator and copropagators 
equations become 

dsQa = +K [V2g„ - N^^q^ (20) 
5.gi = -Aa[v2gt -iVJ^^gt^] (21) 

where V = (l/Lc)V and Aq, = N^/N, and the dimensionless chemical potentials 
are 

M W 

f4 = XI Xij<t)j + X] ^kii^k - (22) 

3=1 k=l 

w 

Af = E4'^^- (23) 

The incompressibility condition V • v = and continuity equations are unchanged 
in dimensionless form 



dtcl>i = -V-ct>ivf (24) 
4Vi = -V-Vivf (25) 
as are the constitutive equations 

&i = Goi^^KT + Ko^Oi^^ " /0(V • Vt)5 " ai/f^, (26) 

if we define the dimensionless shear moduli as Goi = Goi/ac, the dimensionless bulk 
moduli as Kqi = Koi/ac, and the dimensionless relaxation times as Ta = Ta/tc- 
Following convention, we scale each term in the momentum balance equation by 
the characteristic viscous force fc = rjcVc/L'^, where r]c = Cctc, which gives 

= Ca-^ Vtt + Vp-V ■& -f^ -f'f' (27) 

where Ca = is the Capillary number. The equations for the relative velocities 
and wall forces may be expressed 

rcfwi = +afV ■a-4>i (Ca-i V/if + Vp) - if (28) 
if = -aV- V . ^ + ^ . (^Ca-i VA,^ + w) + TC/ (v| - vt) (29) 

where F = is a dimensionless friction factor and = niin(^oi) is the small- 
est monomer friction coefficient. For convenience, we will omit the hats from the 
dimensionless quantities in the subsequent discussion. 



4 Numerical Implementation 

The simultaneous solution of the thermodynamic, hydrodynamic, and viscoelas- 
tic equations requires a fast and efficient implementation to keep the problem nu- 
merically tractable. We propose a multi-step strategy with pseudo-spectral spatial 
discretization and a semi-implicit time discretization which can resolve sharp inter- 
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faces with a minimal number of grid points, and lends itself readily to parallelization 
using freely distributed fftw-mpi fast fourier transform routines. The method is com- 
posed of the following major steps: 

(1) Solve the nf{r) = local thermodynamic equilibrium conditions eq.(lO) using 
an iterative fixed point method and a pseudo-spectral operator splitting scheme to 
obtain the mean field chemical potentials ^'^{r) and /u'^(r). 

(2) Solve the coupled multi-fluid, constitutive equation system using an iterative 
fixed point method which employs semi-implicit time discretization, pseudo-spectral 
spatial derivatives, and the Chorin-Temam projection method [27,26], to obtain ve- 
locity fields v(r) and Wj(r) that simultaneously satisfy momentum balance eq.(15), 
the divergence free condition v = 0, as well as no-flow and no-slip boundary condi- 
tions [20]. 

(3) Use the updated velocities to transport the volume fraction fields ^i(r) and 
in a flux conserving manner. 



4-1 Thermodynamic Iteration 



The chemical potential fields may be obtained by solving local the thermody- 
namic equilibrium conditions fi'^ (r) = which implicitly specify the conjugate fields 
uJj{r) as a function of the volume fractions (pi{v). These conditions are highly non- 
linear, and are satisfied by an iterative procedure which computes the propagators 
and adjusts the conjugate fields to reduce the error. 

The propagator diffusion equations are numerically integrated using a well known 
pseudo-spectral operating splitting technique [18] 



qa{r, s + As)= e-^«^"^V2_;F-i 



Qair, 



(30) 



qi{r,s- As) = e-^«^«^V2_^-i 



-naAaAs/2 t 



(31) 



where J- and J-~^ represent forward and inverse Fourier transform operations, and 
the effective conjugate potential at index s is r2Q(r,s) = ^f"^ ^i(r)7j(s)- The par- 
tition function for each copolymer chain corresponds to the volume average Qa = 
y J ^0(1") 1) and the auxilliary volume fraction fields ^i(r) = ^ dsqa{r, s)qa{r, s)^i{s) 
are computed by quadrature over the polymerization index s. 

The residual errors £i{r) = 0i(r) — 0j(r) are reduced by employing a hill-climbing 
technique with line minimization as described in [19]. A trial step is taken along 
the gradient of the energy surface u>*{r) = (r) -|- eei{r), where e is some small 
adjustable parameter, and a second set of errors e*(r) is calculated using the up- 
dated values. These two points are fit to a parabola to estimate the step size 
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£11— £12 



re 



which minimizes the error in the current search direction, where 
i(r) and £12 = y f drei{r)e*{r). This value is then used to ad- 



vance the fields toward a minimal error solution 



(r) + {a-e)eUr) (32) 

Once the residual error has converged to within some acceptable tolerance, the 
chemical potential fields may be calculated using eqs. (22) and (23). 
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4-2 Hydrodynamic Iteration 



The fields (T, v, Wjjf^ and p are all coupled by the hydrodynamic equations of 
motion and must be obtained simultaneously using an iterative procedure, which 
begins by calculating the updated clastic stresses using a semi-implicit scheme in 
which the nonlinear terms are treated explicitly. 



(33) 



The superscript n counts iteration of the outer, dynamic loop and the superscript 

j counts the iterations of the inner, psuedo-cquilbrium loop. The the linear terms 
are treated implicitly and the upper convected derivatives are computed psuedo- 
spectrally in Fourier space 

— J^~^i (vT • kcTj — kvT ai — a-i ■ vrk) 



rJ - 



as is the shear stress tensor Kj, 



1 + At/Ta 

^-^i [kvT + VTk] - (2/(i)J^- 



ik 



rameter At is the time step chosen for the outer loop. The divergence of the total vis 



(34) 



6. The pa- 



M ~j 



coelastic stress is also computed psueodospectrally, V-<t^ = ^ 
The relative velocity fields are computed using the estimated viscoelastic forces and 
the pressures produced in the previous iteration, 



w; 



J+l _ ^ 



Ca-^V/ip + Vp' 



(35) 



If the walls are externally driven, the motion of each solid object is specified para- 
metrically as a function of time, and the rigid body motions are converted to a set 
of wall velocity fields v^'"(r) which are used to calculate the forces imposed on the 
fluid by the walls. 



a. 



'k 



(36) 



The momentum balance condition is satisfied by seeking a steady state solution for 



over the pseudo-time variable h. 



(37) 



4-3 Chorin-Temam Projection 



The pressure and velocity fields arc obtained simultaneously using the Chorin- 
Temam projection method [27,26]. Gathering the unknown quantities on the left 
hand side of the equation gives 

V* = + /i (v • <7^ - Ca-^VTr" + f^'"^ + f^'^) (38) 

where v* = v-'"'"'^ -|- /tVp'"'"^ may be viewed as the Helmholtz decomposition of a 
single compressible velocity field. Taking the divergence of v* and employing the 
incompressibility condition V-v-'"''^ = produces a Poisson equation for the pressure 
field V^p^"*"^ = V • v*//i which may be solved in Fourier space. 

= -T-^ (^k • v7/ifc2) (39) 
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The divergence-free portion of the velocity field may then be obtained from 

^j+i = V* - hVp>-^^ (40) 



4.4 Volume Fraction Transport 



The preceding sequence repeats until the residual errors, Sy = max || — || 
and £p = max|V-v|, are within acceptable tolerances. Once the hydrodynamic 



iteration has converged, the velocities fields for each component v! 



<f>,n+l 



+ 



i,n+l 



are used to transport the volume fraction fields in a flux conserving manner. 



(41) 
(42) 



4-5 Wall Field Regularization 



Unbounded chemical potentials would be required to produce regions that are 
completely free of a given component. Therefore, in a manner analogous to other 
flow penalization techniques [16], we employ slightly porous walls that do not occupy 
the entire volume at a given location. Solid walls with no slip and no flow boundary 
conditions are recovered in the limit where the porosity is taken to zero as discussed 
in [20]. For the simulations presented in this paper, we employ a porosity of # = 0.01. 

Additionally, pseudo-spectral Fourier methods can produce unwanted Gibbs phe- 
nomena if the system contains overly sharp interfaces. To prevent this, we apply a 
Gaussian filter of the form /(r) = exp(— r^/2a) to the wall fields ipj{r) to smooth 
the transition from solid regions to liquid regions in the channel. The characteristic 
width of the gaussian filter, a, employed in our simulations is typically on the order 
of 0.2Rg, and the filter is applied by convolving the two functions in Fourier space, 

smoothed ~ •^"^(V'j/)- The smallest value of alpha which may be employed is 
dictated by the spatial resolution of the simulation, as the wall transitions must be 
resolved by several grid-points. 



5 Numerical Tests and Validation 

We validate the technique by presenting numerical experiments designed to test 
the capabilities of each model component as well as grid refinement studies which 
demonstrate the accuracy and stability of the numerical implementation. 

5.1 Thermodynamic Properties 

We validate the thermodynamic components by examining phase separating sys- 
tems where hydrodynamic and viscoelastic effects play a minimal role. In particular, 
we consider the phase separation of diblock copolymer melts in two and three di- 
mensions. This system is of great interest for patterning applications as it exhibits 
geometrically regular repeating structures with nanometer length scales resulting 
from the competition between enthalpic and entropic forces. The structure and 
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Fig. 2. Equilibrium diblock copolymer meso-phases are obtained in the long time 
limit, (a) Lamellae, /a = 0.5, xN = 18 (b) Gyroid, /a = 0.37, xN = 18 {QiaM 
space group) (c) Hexagonally packed cylinders, /a = 0.3, x^ = 21 {QimSm space 
group) (d) Close packed spheres, /a = 0.23, x^ = 0.20 

phase boundaries of these state have been examined in great detail using equi- 
librium SCFT techniques, and while HSCFT is a non-equilibrium formulation, it 
should produce the same equilibrium morphologies in the limit of long simulation 
times. 

The equilibrium state exhibited by a particular system is determined by the vol- 
ume fraction of the components /j and the degree of incompatibility x^ between 
them. In fig. 2, we demonstrate that the HSCFT method produces the anticipated 
equilibrium morphologies in three dimensions including lamellae, gyroids, hexago- 
nally packed cylinders, and close packed spheres. Comparison with the phase bound- 
aries calculated in [25] confirms that the morphologies obtained are consistent with 
their coordinates in phase space. 

We examine the dynamic evolution of this same system to test the numerical 
convergence of the technique under grid-size refinement. Consider a symmetric block 
copolymer melt where both monomer species occupy an equal volume fraction, 
fA = fB = 0.5, and the Flory interaction strength, NxAB = 20, places the system in 
the intermediate segregation regime. The melt is initialized to a homogeneous state 
with small, random fluctuation imposed at all frequencies to break the symmetry. 
Fig. 3 illustrates the time evolution of this system after it has been rapidly quenched 
below its order disorder transition temperature. Initially, high frequency modes 
rapidly decay, and density fluctuations grow exponentially, with the greatest growth 
occurring for structures at a critical wave number. The fluctuations saturate, leaving 
the melt in a highly defect-filled meso-phase separated state in which the feature 
size grows over time, until the connectedness of the diblock copolymers prevents 
further coarsening. Morphological defects may persist for long times, correspond to 
kinetically trapped states. 

We use the same pattern of random variations about the homogeneous state 
to initialize simulations of 32x32, 64x64, and 128x128 grid-points in a periodic 
system of fixed dimension, 16 x IQRg. When we superimpose the volume fraction 
contours at all three simulations, we find that the higher resolution simulations 
converge to a single solution as the mesh is refined. However, the lowest resolution 
predicts a distinctly different defect pattern after the fluctuations saturate. We 
conclude, therefore, that at intermediate segregation strengths, two gridpoints per 
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Fig. 3. Convergence of diblock copolymer meso-phase separation under grid refine- 
ment, [left]: /a = 0.5, = 20 at t = 0, 1, 3, 10, and 80. [right]: /a = 0.3, = 20 
at t = 0, 1, 4, 10 and 120. Snapshots where taken at three resolutions (a) 32x32 (red), 
(b) 64x64 (green), (c) 128x128 (blue), (d) Superposition of meso-phase boundaries 
at all three resolutions at f = 4. 

Rg is insufficient to fully resolve the mesoscale structures, but four or eight grid 
points is adequate. 

5.2 Hydrodynamic Properties 

We validate the hydrodynamic component by simulating low Reynolds number, 
pressure driven and wall driven flows in the absence of thermodynamic or viscoelas- 
tic effects. Three classic examples of Stokes flow found in many texts [24] include 
flow past a fixed sphere, flow induced by a moving sphere, and the flow generated 
in a journal bearing due to the motion of two co-rotating eccentric cylinders. The 
streamlines produced in the simulation of each case is consistent with the known 
analytical predictions as illustrated in fig. 4. The flow past a fixed sphere (a) exhibits 
the expected azimuthal symmetry and time reversibility. The streamlines diverge as 
they flow past the fixed sphere (a) and converge in the presence of the moving sphere 
in fig. 4(b). The third case (c) corresponds to the counter clockwise rotation of an 
infinite cylinder in a fixed cylindrical cavity, which produces a counter clockwise 
rotation in the region closest to the rotating cylinder and a reversed, clockwise flow 
in the regions furthest from it. This result is also consistent with analytic solutions 
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Fig. 4. Streamlines of wall driven and pressure driven viscous flows at low Reynolds 
numbers with periodic boundary conditions, (a) Flow past a fixed sphere (b) Flow 
induced by a moving sphere (c) Flow reversal in a journal bearing 

of the Stokes equations [24] as well as with experimental observations [23] . 

Wc examine the pressure driven flow past a flxed sphere in greater detail to 
test the behavior of the hydrodynamic implementation under grid refinement. The 
velocity field is generated by an external force, = 0.2 which drives a melt of 
homopolymer A from top to bottom in the channel. The system also contains a 
transverse band of incompatible homopolymer B used to track the fluid motion. The 
Flory incompatibility is set at Nxab = 8, and the simulation domain is periodic 
with dimension 10 x 2QRg. The fixed sphere occupies 40% of the channel width and 
interacts neutrally with both materials. The viscosity is r]s = 1.0, and viscoelastic 
effects are suppressed, Goi = Kqi = 0. 

The interaction of the transverse homopolymer band with the fixed spherical ob- 
stacle produces a sequence of events in which the homopolymer band wraps around 
the sphere and is stretched until it ruptures, leaving a layer of B in contact with 
the obstacle and a drop of homopolymer B which coalesces in the wake of the ob- 
stacle, as shown in fig. 5. We repeat the simulation at grids resolutions of 32x64, 
48x96, and 64x128, and compare the volume fraction contours produced by each 
in fig. 5a-d. The contours of homopolymer A are superimposed and are found to 
be nearly identical, indicating that the method converges to a single solution as 
the simulation mesh is refined. We also characterize the fiow rate by measuring the 
maximum Weissenberg number Wi = rfmax/-^ where fmax is the maximum mean 
field velocity, and L is the width of the narrowest portion of the channel. In all three 
cases, we obtain the same value, Wi = 0.32. In fig. 5(e) we analyze this system in 
three dimensions illustrating that the method produces results consistent with the 
two dimensional simulations. 

5.3 Viscoelastic Properties 

We validate the viscoelastic components of the model by examining the viscoelas- 
tic phase separation of two component systems. "Viscoelastic phase separation" is 
the term given to the de-mixing process in two fluids systems which exhibit a large 
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Fig. 5. Convergence of the viscous, low Reynolds number flow of a band of ho- 
mopolymer B past a circular cylinder in a matrix of homopolymer A under grid 
refinement. Snapshots are taken at t = 0, 2, 4, 6, 8, and 10 at three resolutions, (a) 
32x64 (red) (b) 48x96 (green) (c) 64x128 (blue), (d) Superposition of homopolymer 
A interface at all three resolutions at t=6. (e) The same system in three dimensions. 

contrast in their dynamic properties (dynamic asymmetry.) In such a system, phase 
separation occurs in a manner which differs from the usual spinodal decomposition 
process. Tanaka observed this behavior in polymer solutions [22] and modified the 
two fluid model to account for them [2] . His model has subsequently been applied 
to viscoelastic phase separation of polymer blends and diblock copolymer melts as 
well [1,21]. 

In fig. 6, we illustrate the time dependent behavior of six systems which ex- 
hibit phase separation behavior. Cases (a) and (f) are included for reference, where 
case (a) illustrates normal phase separation in a critical diblock copolymer melt, 
/a = 0.5, and case (f) demonstrates phase separation in an off-critical melt Ja = 
0.30. Cases (b)-(d) demonstrate viscoelastic phase separation behavior of systems 
in which a higher bulk modulus material Ka = 5, = 2 (white) has been mixed 
with a softer material Kb = 0, = 1 (black) both of which demonstrate nonzero 
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Fig. 6. Viscoelastic phase separation of two homopolymer components, (a) Normal 
AB diblock, Ja = 0.5 (b) Glassy AB diblock, /a = 0.5 (c) Glassy A+B blend, 
/a = 0.5 (d) Glassy A+B blend, /a = 0.3 (e) Glassy AB diblock, /a = 0.3 (f) 
Normal AB diblock, /a = 0.3 



shear moduli Ga = Gb = 0- All four systems exhibit the expected characteristic 
behavior beginning with an incubation period in which the hard phase forms a vis- 
coelastic network structure, followed by nucleation of the softer phase, and eventual 
network break-up. Case (b) is a critical /a = 0.5 glassy /elastic diblock and case (c) 
is a critical blend. Case (d) is an off-critical ,fA = 0.3, glassy/elastic blend and (e) 
is an off-critical diblock melt. As expected, after the breakup of the network, each 
material resumes its phase separation processes in which the blends macro-phase 
separate and the diblocks do not. Cases (d) and (e) demonstrate "phase inversion" 
wherein the soft material forms drops in a hard matrix in the early stages and their 
roles are reversed in the late stages. Each of these results is consistent with previous 
investigations [1,21] 
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6 Conclusions 



To summarize, we have introduced a method caUed Hydrodynamic Self-Consistent 
Field Theory (HSCFT) which extends the capabilities of SCFT to the hydrodynamic 
regime in which bulk material transport and viscoelastic effects play a significant 
role. We introduced a semi- implicit, iterative numerical scheme for the solution 
of these equations which makes liberal use of pseudo-spectral fast-Fourier tech- 
niques, enabling practical simulation of sizable systems. We validated the thermo- 
dynamic component by reproducing the correct cquilibriTim mcso-phases found in 
three dimensional diblock copolymer melts. We validated the hydrodynamic model 
by demonstrating that it produces low Reynolds number flows that are consistent 
with analytic solutions of the Stokes flow equations. The viscoelastic constitutive 
equations were validated by simulating the expected viscoelastic phase separation 
dynamics of glassy /clastic AB diblock copolymer melts and glassy/elastic A-|-B 
blends. We also performed a series of resolution studies on systems with and with- 
out bulk flow in which we demonstrated that the combined system is stable and 
converges smoothly to a single solution under grid-size refinement. We believe that 
the coupled interaction of this system of equations provides a powerful and flexible 
tool for studying the behavior of complex fluids in hydrodynamic flows. 

A Derivation of the Thermodynamic Model 

In this section, we detail the derivation of a non-equilibrium, self consistent 
field theoretic model general enough to represent an arbitrary blend of C distinct 
multiblock copolymers with varying molecular weights. The derivation follows the 
standard procedure outlined in [13], which consists of the following steps: 

(1) Construct a particulate, mesoscale model of the essential physics. 

(2) Convert the model to field theoretic form. 

(3) Approximate the partition function and obtain the thermodynamic forces. 
A.l Mesoscale Model Construction 

The most general system we arc interested in describing consists of C copolymer 
species and W solid wall materials, which may be visualized as a dense array of 
mesoscale particles or "monomers" of equal volume vq. The forces between the 
monomers consist of a hard core repulsion, covalent bonds between monomers on 
the same polymer, and Van der Waals forces between non-bonded monomers. 

The hard core potential is modeled implicitly by enforcing the incompressibility of 
the copolymer melt such that the total number density remains constant X^^i Pi = 
po = 1/vq. The sum runs over all distinguishable monomers M and the number 
density operator for monomer i is Pi{r) = Yl^=i [q " ~ Raj(s)) 7j(s)], where 

the quantity Ua counts the number of copies of copolymer type a and Na is its 
polymerization index. 

Covalent bonds are modeled by treating each group of wall particles as a single 
rigidly rotating and translating object and each copolymer as a string of monomers 
on a parameterized space curve Rai(s) where a G C indicates the copolymer type 
and i E na- Using the standard Gaussian thread approximation, the unperturbed 
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polymer chains are represented as random walks with a "stretching" free energy 
given by 



kT 



C Ua 

a=l i=\ 90c 
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Jo 



ds 



dHaiis) 



ds 



(A.l) 
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where Rga is the unperturbed radius of gyration of copolymer species a. 

Van der Waals interactions between non-bonded particles are represented by the 
interaction potentials 

M M 

i=i j=i 

M W 
i=i j=i 

where U^^ is the total energy of monomer-monomer interactions, [7^^ is the total 
energy of monomer-wall interactions, and the volume fraction operator is defined as 
0i(r) = pi{r)/po. Similarly, V'j(r) represents the volume fraction operator of solid 
material j at r. The strength of the effective repulsion between dissimilar units is 
given by the Flory parameters x a-^d 

Wc complete the mesoscale model by constructing a partition function over all 
physically realizable configurations of the system in the canonical ensemble 



kT 



Z = ZoJ Y[V[Ilai]V[r]d{^j{r) - ^j{r))eM-U/kT) 

oi,i,j 



(A.4) 



where the total energy of a given configuration is ?7[{Rai}] = Uq + U^cf) + U^^. 
The ensemble sum runs over all possible space curves for each copolymer with 
delta functions which select the configurations that are consistent with the volume 
fraction fields ^^'i(r). 

A. 2 Conversion to Field Theoretic Form 



The mesoscale model is converted to field theoretic form using a formal procedure 
which replaces the particle degrees of freedom with a set of continuous volume 
fraction and chemical potential fields. Using the delta functions, the volume fraction 
operators 0(r) are replaced with their equivalent field values ^(r) in the interaction 
terms 
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(A.5) 



and the conjugate chemical fields iOi{r) are introduced by employing a Fourier rep- 
resentation for the delta functions. 



P[r]5(</.,(r)-,^,(r)) 



exp 
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Vo 



drujj{r)[(i)j{r) - (t>j{r)] 



(A.6) 



Because each instance of a copolymer of type a is indistinguishable from others 
of the same species, the part of Z which depends explicitly on R^j factors into a 
product of Ua identical single chain partition functions when we insert the definition 
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of the volume fraction operators <j) 



Qa — Qo 
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ds 





(A.7) 



where we have defined fia(r, s) = Yli^i <^i{''^)li{^) ■ 

This integral over all paths may be converted into a volume integral over a 
propagator in a manner analogous to that employed by Feynman and Kac [14] in 
their path integral formulation of quantum mechanics. When expressed in this form, 
Qa = y J drqa{r,s)qlt{r,s), and the propagator and co-propagator are evaluated 
by numerically solving the differential equations with initial conditions ^^(r, 0) = 1 
and qi{r, 1) = 1. 

dsQair, s) = +Rl^V\a{r, s) - Nana{s)qa{r, s) (A.8) 
dsqiir, s) = -Rl^V^qiir, s) + 7V,0,(s)9t (r, s) (A.9) 

The partition function Z may now be expressed in the purely field theoretic form 
Z = Zq J Yl- VluJi] exp (F[{u;, (/)}]/kT) where the quantity F may be interpreted as 
the free energy associated with a particular configuration of fields uj(r) and 0(r) 
which takes the form 
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M W 
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a 

(A.IO) 



A. 3 Mean Field Approximation 

Since we do not know how to evaluate the partition function, it must be nu- 
merically sampled or analytically approximated. As discussed in [13], the chemical 
potential fields uji are in general complex, which makes direct numerical sampling 
difficult and time consuming. Instead, a mean field approximation of the partition 
function is obtained by performing a stationary point analysis on the integral for 
the conjugate fields Wi(r). The resulting solution is purely real and may be obtained 
from the local thermodynamic equilibrium conditions (LTE) 

JI- = ^[^,{r)-Ur)]=0 (A.ll) 
du;i{r) vq 

The quantitie ^i(r) = — t^oWQ ^^^^ is called the auxiliary monomer volume fraction, 
which may be expressed as an integral over the propagators 

^^(^) = 7T- f dsqa{r,s)ql{r,sH{s) (A.12) 



Qa Jo 

where = naNa{vo/V) is the volume fraction occupied by all copolymers of species 
a in the system. 

The LTE conditions represent a set of highly nonlinear equations which spec- 
ify the mean field conjugate potentials a;j(r) as a function of the volume fraction 
fields (j)i{r). The conjugate fields correspond to local configurations of the copoly- 
mer chains which relax quickly in comparison with the conserved volume fraction 
fields, and as such are be expected to fluctuate close to their mean field values. Once 
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we have solved the LTE conditions, we may compute the non-equilibrium chemical 
potentials fi'^{r) and /i^(r) which correspond to functional derivatives of the free 
energy with respect to the volume fraction fields. 

^'ti^^ = = E xMr) + E ^Mr) - u^iir) (A.13) 
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B Derivation of the Multifluid Model 

In this section, we derive hydrodynamic equations of motion for the transport 
of multiple viscoelastic fluids in the presence of rigid channel walls using Rayleigh's 
variational principle as discussed in [15]. We briefly review the method they em- 
ployed to derive a two-fluid model for polymer blends and apply the same procedure 
to derive its multi-fluid generalization. 



B. 1 Review of Rayleigh 's Variational Principle 

For systems that are not too far out of equilbrium, the theory of irreversible ther- 
modynamics makes the approximation that the generalized velocities in the system 
are linearly related to the generalized forces by ^ = — Ylj ^ij M~ where Xi are the 
generalized systems variables, and are Onsager coefficients. For cases where the 
inverse of the kinetic Onsager coefficient matrix is defined, this relationship may be 
re-expressed as ^ -|- L^-^^ = 0. This equation of motion may be equally well 
expressed in terms of a variational principle in the generalized velocities SR = 
where R = §^ii + \ Y^ij^7j'^i^3- The two terms in the Rayleigh functional 
may interpreted as the total change in free energy F and an energy 
dissipation function W . 



B.2 Multi-fluid Model Derivation 



The total change in free energy for the system under consideration is 



dr 
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(B.l) 



In the absence of chemical reactions, the number density of each component is a 
conserved quantity such that the monomer fields obey 0j = —V • 4>ivf and the 
solid fields obey ipj = —V • V'jV^- Upon substitution of these relationships, the total 
change in free energy may be written in terms of the monomer velocities vf and 
the rigid wall velocities as 
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(B.2) 



To construct the dissipation function W, we note that energy is dissipated due 
to friction between the monomers and the entangled polymer network Ci^(vf — vt)^, 
by friction between the solid material and the network Cfi'vf — v^)^, and by elastic 
deformation of the network cr : Vvj^. Energy is also added to the system by the 
motion of the externally driven rigid walls fj^ • and by body forces acting on the 

bulk of the fluid, if ■ ^r'f . Therefore, we construct a dissipation function of the form 



W_ 
~2 



(B.3) 



The velocity field describing the motion of the entangled polymer network is 
called the tube velocity = Yli ^t^t X^j ^^^i which is a rheological mean 
velocity where each component is weighted by the stress division parameters. The 
tube velocity vr may also be expressed in terms of the mean velocity v and the 
relative velocities Wj = vf — vt using the following relation. 



vr 



E,K-^Owi + E,-K-V'>f + v 



1 - E.K - 4>k) 

Solving a force balance condition on the entangled network Cfi^t 



(B.4) 



Y^j Cj{^t ~ Vt) = produces the stress division parmeters 



Vt) + 
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C//C where the sum of the friction coefficients is ^ = Cf + Cf- The friction 

coefficients take the form (f = Coj(^a/-^ea)</*t as discussed in [15] where Nea is the 
entanglement length of polymer species a. The friction coefficient for flow through 

Cojipj/^ as discussed in [20,16], 



a semi-porous wall is given by Darcy's law, 



where the porosity of the material is $ = YlT = 1-0 — YlY ^i' 

The equations of motion for each component are obtained by appending the 
divergence free condition to the Rayleigh functional R = F + — p(V • v) and 
then extrcmizing it with respect to the velocity fields, = and = 0. The 

full Rayleigh functional (with implied summation over repeated indices) is 
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which produces the following equations of motion, 
= cf>,v4 + cf>,Vp - afv ■ a + C/(vf - vt) - 



0: 



(B.6) 

: ^jV/x| + ^jVp - ap ■ a + C/(v,^ - vt) - f/ (B.7) 

Eq. (B.6) may be solved for the velocity of each fluid relative to the elastic polymer 
network Wj = — vt, and Eq.(B.7) gives the local forces fj'(i') that the walls 
impart upon the fluid when moving at a velocity vj^ (r) . 
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(B.8) 



:V.,V^,^ + V',Vp-4V-cT + C/(v^^-vr) (B.9) 

Summing Eq. (B.6) and Eq. (B.7) over all components, gives the momentum balance 
condition 



= Vvr + Vp - V • cr - - f 



(B.IO) 



in the limit of zero Reynolds number, which together with the divergence free 
condition implicitly specifies v. The total osmotic pressure is defined as Vtt = 
J2f (t)iVfif+J2Y The net external force exerted by the walls is = J2Y 4 

and the total body force acting on the fluid is = J2i^ ■ 
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